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Time-reversible Born-Oppenheimer molecular dynamics 
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We present a time-reversible Born-Oppenheimer molecular dynamics scheme, based on self- 
consistent Hartree-Fock or density functional theory, where both the nuclear and the electronic 
degrees of freedom are propagated in time. We show how a time-reversible adiabatic propagation of 
the electronic degrees of freedom is possible despite the non-linearity and incompleteness of the self- 
consistent field procedure. Time-reversal symmetry excludes a systematic long-term energy drift for 
a microcanonical ensemble and the number of self-consistency cycles can be kept low (often only 2-4 
cycles per nuclear time step) thanks to a good initial guess given by the adiabatic propagation of the 
electronic degrees of freedom. The time-reversible Born-Oppenheimer molecular dynamics scheme 
therefore combines a low computational cost with a physically correct time-reversible representation 
of the dynamics, which preserves a detailed balance between propagation forwards and backwards 
in time. 

PACS numbers: 71.15.Pd,31.15.Ew,31.15.Qg,34.10.-|-x 



Ab initio molecular dynamics based on Hartree-Fock 
or density functional theory [illllllllll[llili 

has become an important tool for simulations of an in- 
creasingly wider range of problems in geology, material 
science, chemistry, and biology. Ab initio molecular dy- 
namics, where the atomic positions move along classical 
trajectories, can be categorized in two major groups: La- 
grangian Car-Parrinello molecular dynamics and Born- 
Oppenheimer molecular dynamics [1, H, Ifl H Hi B U^ 
fill lia | . The tremendous success of Car-Parinello molec- 
ular dynamics, invented two decades ago |3|, is based 
on its low computational cost combined with the con- 
served Lagrangian properties of the dynamics. However, 
unless a Car-Parinello simulation is performed carefully, 
it may yield results different from Born-Oppenheimer 
molecular dynamics [3, M HM HI; HB UM- ^"^ Born- 
Oppenheimer molecular dynamics the atomic positions 
are propagated by forces that are calculated at the self- 
consistent electronic ground state for each instantaneous 
arrangement of the ions. Born-Oppenheimer molecu- 
lar dynamics is computationally expensive compared to 
Car-Parinello dynamics because of the requirement to 
reach a self-consistent field (SCF) solution in each time- 
step. However, the number of SCF cycles and thus 
the computational cost can be strongly reduced by us- 
ing an initial guess for the electronic degrees of freedom 
p{tn+i) (here represented by the electron density), which 
is given by an extrapolation from previous time steps 
d, H, y, US Ei Ii3 • The electronic extrapolation scheme, 
combined with the SCF procedure, can be seen as an adi- 
abatic propagation of the electronic degrees of freedom. 



where 



p(t„+i) = SCF[/5(i„),P(in-l) 



(1) 
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Unfortunately, this approach has a fundamental prob- 
lem. Because of the non-linear and irreversible SCF pro- 
cedure, which in practice never is complete, the time- 
reversal symmetry of the electronic propagation is broken. 
This problem does not occur in Lagrangian Car-Parinello 
molecular dynamics |2Clj| , where both the nuclear and the 
electronic degrees of freedom can be propagated with 
time-reversible Verlet integrators |2l| . The main purpose 
of this letter is to show how an effective time-reversible 
propagation of the electronic degrees of freedom is pos- 
sible in Born-Oppenheimer molecular dynamics, despite 
an irreversible and approximate SCF procedure. 

Computational schemes using time-reversible integra- 
tors give an energy meandering around a constant value 
that does not drift with time. This stability follows from 
the time-reversal symmetry, which excludes a steady in- 
crease or decrease of the energy (over the Poincare' time) 
for a periodic motion. Because of the broken time- 
reversal symmetry in the adiabatic propagation of the 
electron density in Eq. (Q), a small but systematic en- 
ergy drift occurs in the evolution of a microcanonical en- 
semble, with an accumulating phase space error. The er- 
ror can be systematically reduced by improving the SCF 
convergence or be removed completely by using an ini- 
tial guess in the SCF procedure that is independent of 
previous time steps 8j _9J . Both these remedies are com- 
putationally expensive and using an initial guess that is 
independent from previous time steps often require a sig- 
nificantly increased number of SCF cycles. 

The basic principles for the time-reversible lossless in- 
tegration of the electronic degrees of freedom is described 
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FIG. 1: The principle behind the time-reversible lossless fil- 
ter process for propagation of the electronic density. In the 
forward filter process (upper part), Pn+i — —pn-i + U{pn) 
and in the lossless backward reconstruction (lower part), 
—pn-i ~ Pn+1 — U{p„). The dual propagation with the 
auxiliary density p„ allows a perfect reconstruction of the 
self-consistent Born-Oppenheimer density p„, despite an irre- 
versible, incomplete and approximate SCF procedure. 



here in terms of the propagation of the electron den- 
sity p{t). However, time- reversible Born-Oppenheimer 
trajectories can be constructed by replacing the den- 
sity by other parameters governing the electronic degrees 
of freedom, such as the effective single-particle poten- 
tial, Hamiltonian, density matrix, or wavefunctions. Our 
approach is therefore general and applicable to a large 
number of electronic structure schemes based on self- 
consistent Hartree-Fock or density functional theory. 

A reversible propagation of the electron density p{t) in 
finite time steps of 5t (i„ = io + nSt) can be constructed 
from a lossless filter process analogous to, for example, 
the lossless wavelet transform used in data compression 
p2|. The principle of the process, which is the key re- 
sult in this letter, is shown in the upper part of Fig. ^ 
The scheme is divided into two channels; the upper chan- 
nel, with an approximate auxiliary density (denoted by 
a tilde), Pn+i = p{tn+i), is used as an initial guess for 
the Born-Oppenheimer density, Pn+i, in the lower chan- 
nel. The Born-Oppenheimer density is given through the 
non-linear, in practice incomplete, and numerically lossy 
SCF procedure, 



Pn+l = SCF[p„+i]. 



(2) 



The function U{pn), which is allowed to be numerically 
approximate and without a unique inverse, is an update 
filter for the propagation of the auxiliary density, 
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U{pn 
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(3) 



It is easy to see that the filter process is perfectly lossless 
and reversible by running the process backwards in time 
with the © sign replaced by a sign, as shown in the 
lower part of Fig. ^ The scheme is therefore a bijective 



map which allows perfect reconstruction of p„ backwards 
in time, despite the fact that the SCF procedure by itself 
is an irreversible lossy transform. 

The auxiliary density Pn+i in Eq. ^ will be close to 
the self-consistent Born-Oppenheimer density Pn+i if the 
lossless filter process in Fig. ^ approximates the time- 
reversible adiabatic evolution of the density on the Born- 
Oppenheimer potential energy surface. This reduces the 
number of SCF cycles necessary to reach the new self- 
consistent density. One way to achieve this is to construct 
the update filter U(pn) in Eq. ^ from the time-reversible 
Verlet integrator [2j such that 



Pn+l — [^Pn 



St' 



Pn\ - Pn-l- 



(4) 



This integrator fulfills time-reversal symmetry since it 
remains the same if we switch the sign of 5t and thus 
interchange Pn-i with Pn-i-i- Note that perfect lossless 
reconstruction is a necessary, but not a sufficient, con- 
dition for time-reversal symmetry. The simplest time- 
reversible approximation of the second order time deriva- 
tive, Pn = d'^pn/dt', is to set it equal to zero. In this case 
the lossless propagation of the auxiliary density is 



Pn+l — "^Pn — Pn-l- 



(5) 



This surprisingly simple difference approximation fulfills 
time-reversal symmetry, allows a perfect reconstruction 
backwards in time, and approximates the propagation of 
the density with a modified linear interpolation from two 
previous time steps. It also avoids an unstable exponen- 
tial error growth since the characteristic equation (with 
Pn replaced by p„) has no roots outside the unit circle. 
If the auxiliary density Pn-i in Eq. |(SJ) is replaced by 
the self-consistent Born-Oppenheimer density p„_i the 
propagation scheme is identical to a linear interpolation. 
Below we demonstrate how such a modification affects 
phase-space conservation and the energy drift. 

A quite general scheme for constructing efficient time- 
reversible integrators for the approximate adiabatic prop- 
agation of the electronic degrees of freedom is given by a 
least square fit of the ansatz 
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pi-t), 



(6) 



m=0 



to Born-Oppenheimer densities p(t„) at N successive 
time steps. A similar least square approximation, but 
without the constraint of time-reversibility, was recently 
proposed by Pulay and Fogarasi for the extrapolation of 
the single-particle Hamiltonian in a highly efficient Fock 
matrix dynamics (FMD) method |a, Q- Often only 2-3 
SCF cycles are necessary in their scheme, but because of 
the broken time-reversal symmetry a small but system- 
atic energy drift occurs. 

By choosing different numbers of fitted values N we 
can calculate the expansion coefficients am in Eq. @ and 
express them in terms of previous Born-Oppenheimer 
densities. For example, a least square fit using the ansatz 
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FIG. 2; The fluctuations in the total energy as a function 
of time for a F2 molecule using Hartree-Fock theory with a 
Gaussian basis set (RHF/3-21G). The time-reversible prop- 
agation based on Eq. ^, with the density replaced by the 
effective single-particle Hamiltonian, is compared to the en- 
ergy drift using the corresponding linear interpolation from 
previous time steps. The (4-2) Fock matrix dynamics (FDM) 
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FIG. 3: Detail of the phase space during 500 fs for one of the 
atoms in the F2 molecule (Hartree-Fock theory with a Gaus- 
sian basis set RHF/3-21G, 2 SCF/step). The time-reversible 
propagation preserves the phase space whereas the non-time- 
reversible dynamics has a small but noticeable drift. 

in Eq. lO, with M — I for 6 previous densities, leads to 
the time-reversible approximation 

Pn+l = T:7[30(p„+p„_4)-3(p„_i+p„_3)-28/9„_2]-Pn-5- 

(7) 
Note that the least square ansatz in Eq. ^ also can be 
extended to negative time exponents that are even. 

To demonstrate the time-reversible lossless Born- 
Oppenheimer molecular dynamics we have implemented 
the scheme in MondoSCF J23 , a suite of programs using 
Gaussian basis sets for electronic structure calculations 
based on self-consistent Hartree-Fock or density func- 
tional theory. The density in the time-reversible prop- 
agation in Eqs. Q and IjHJl has been replaced by the ef- 



fective single-particle Hamiltonian. The number of SCF 
cycles has been measured in the number of constructions 
of Hamiltonians, which is the most time consuming step. 

Figure|2]shows the total energy as a function of time for 
an F2 molecule. The modification of the time-reversible 
linear propagation in Eq. @ to a linear interpolation 
leads to a systematic drift in the energy. Using more 
SCF cycles per time step reduces the energy drift, but it 
never really disappears. As a comparison we also show 
the (4-2) Fock matrix dynamics (FMD) scheme by Pu- 
lay and Fogarasi iS] , based on a second order polynomial 
least square fit using four previous data points. This 
scheme gives, in principle, a more accurate extrapolation, 
which is noticed in a smaller amplitude of the energy os- 
cillations. However, because of the broken time-reversal 
symmetry there is a systematic drift in the energy. For 
the time-reversible linear integrator in Eq. (0), using only 
2 SCF iterations per time step, any energy drift was less 
than 10~® Hartree/ps. The phase space is also conserved 
with the time-reversible integration, as illustrated in Fig. 
P|). Thus, even with an incomplete SCF convergence 
time-reversibility is correctly preserved. 

Though unproven, all time-reversible propagators we 
have found so far have characteristic equations with all 
their roots on the unit circle. Because of the perfect loss- 
less reconstruction, any error that occurs in the calcu- 
lations will propagate throughout the simulation. This 
leads to a random noise that increases with time. Be- 
cause of this noise the auxiliary density pn slowly moves 
away from the self-consistent solution. This means that 
it is not possible to simply take long time steps using only 
1 SCF cycle per time step. An increased number of SCF 
cycles (or shorter time steps) is in general necessary to 
reach a sufficiently accurate Born-Oppenheimer density 
for longer simulation times. Figure 0] shows the fluctua- 
tions in the total energy for a C2F4 molecule during 1 ps 
of simulation time at a temperature T « 500K. The ex- 
trapolation scheme in Eq. Q was used and 3 SCF cycles 
per time step was applied using Pulay's direct inversion 
in the iterative subspace (DIIS) algorithm to accelerate 
the convergence [24ll25| . 

The ability of a perfect reconstruction of the dynamical 
data backwards in time is kept also using approximate 
arithmetics, thanks to the lossless bijective filter process 
illustrated in Fig. ^ The lossless property is therefore 
kept also when small elements below some drop tolerance 
are set to zero in the SCF optimization of p„ and in the 
update function C/(p„). This is potentially important 
in the study of very large systems using, for example, 
linear scaling electronic structure methods |2g| that in 
general require a lower numerical accuracy compared to 
conventional schemes. 

The main result in this letter is that we have shown 
that an effective time-reversible propagation of the elec- 
tronic degrees of freedom is possible in self-consistent 
Born-Oppenheimer molecular dynamics, despite the irre- 
versible, non-linear, and in practice always approximate 
SCF procedure. This ability may also open the way for 
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FIG. 4: The energy fluctuations around the average energy 
as a function of time for a C2F4 molecule at a temperature 
T ~ 500K (Hartree-Fock theory with a Gaussian basis set 
RHF/3-21G). Three SCF cycles per time step of length St = 
0.25 fs were used. 



the development of higher order symplecitc integrators 



[23,12^123,133,133 in ab initio molecular dynamics, 

In summary, we have presented a scheme for time- 
reversible Born-Oppenheimer molecular dynamics that 
combines a low computational cost with a physically 
correct time-reversible representation of the dynamics. 
The principle is based on a lossless filter integration of 
the electronic degrees of freedom, where a dual prop- 
agation with an auxiliary density allows a perfect re- 
construction backwards in time of the self-consistent 
Born-Oppenheimer density. Time-reversal symmetry re- 
moves a systematic long term energy drift for a micro- 
canonical ensemble. This energy stabilization, combined 
with the adiabatic propagation of an approximate Born- 
Oppenheimer electron density, strongly reduces the num- 
ber of necessary self-consistency cycles, while maintain- 
ing an accurate physical description. 

This work was performed under the auspices of the 
US Department of Energy, Office of Science and LANL 
LDRD/ER program. A.N. acknowledge support form the 
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